%% master_counterfactuals.m

%clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Date: February 2019
%
% Purpose: This file runs the counterfactual exercises after the
% main_baseline.m has been run. Change the shocks to experiment with the
% counterfactuals.
% 
% Based on one file: Step_4_Algorithm_counterfactual.m, which uses new
% baseline input data, MAT/Step_4_setup_`alphaT''rhoT'_`shockT'.mat
% as initial level variables, and then shocks the system and computes hat
% variables for wages, prices, exports, etc.
%
% NEW: for CES setup. Have streamline file that has to be called up and
%                     various variables that have to be used.
% Output: MAT/Step_4_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat,
%         MAT/Step_5_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%         MAT/Elasticity_output_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma lambda eta ...
    tol maxit D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% alphaT: 
%
% 1 = PWT labor share for all countries and Labor share for French firms 
%     are at the COUNTRY level.
%
% 2 = PWT labor share for all countries and WIOD scaled labor share for 
%     French firms. Labor share for French firms are at the FIRM level.
%
% 3 = WIOD labor share for all SECTORS with the labor share at the SECTOR 
%     level for French firms.
%
% 4 = WIOD labor share for all SECTORS with the labor share at the FIRM 
%     level for French firms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Loading in the previous data to get some useful variables

fname = strcat('MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat');
load(fname,'FI_J', 'firmsecD_sorted', 'start_sorted', 'alphaT');
load('MAT/WIOD.mat','countrysecD','france','J','N','start','sum_by_country_dummy');

%% SET PARAMETERS VALUES

% Frisch elasticity for GHH preferences

%varphi=3;

% share of labour in total final consumption (need to make assumption or collect data)

sL_n = ones(1,N)*0.8; 

% share of profits in total final consumption (should have this data)

sPi_n= ones(1,N)*0.15; 

% share of trade deficit in total final consumption (should have this data)

sD_n= ones(1,N)*0.05; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rhoT:
%
% Value for rho. Will apply to J number of sectors to create a vector of
% homogeneous rhos across sectors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% rhoT = '5';
% 
% rho = 5*ones(J,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Value for sigma
%
% This is the elasticity of substitution across sector consumption. Given
% we do not adjust this anymore after running C-D production we set to the
% value we think is most reasonable: 1.5. There is therefore no longer a
% "type" variable.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sigmaT = '1.5';
% sigma = 1.5*ones(J,1);
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lambdaT:
%
% Value for lambda. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% lambdaT = '1.000001';
% 
% lambda = 1.000001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etaT:
%
% Value for eta. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% etaT = '1.000001';
% 
% eta = 1.000001;

%% Baseline: no shocks

% Pi_hat: Change in profits per country
    
Pi_hat_n=ones(1,N);
    
% D_hat: Change in trade deficit
    
D_hat_n = ones(1,N);
    
% Taste shock - for France only
    
Xi_hat_mnjf=ones(FI_J,N);
    
%Taste shock - for sectors only (countries except France)
    
Xi_hat_mnj= ones(N*J,N);
    
% Productivity shock, a is firm specific - France only
    
a_hat_f= ones(FI_J,1);
    
% Productivity shock, a is for sectors - countries except France
    
a_hat = ones(J*N,1);

for qq = 2:2:4
    %1:4
    
%% Feeding in SHOCKS
    
    if qq == 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 1:
    %
    % A +10% productivity shock to all firms for Germany. Germany is the 10th
    % country in the WIOD.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'DEUprod_p10'
    
    a_hat = ones(J*N,1);
    a_hat(start(10):start(10+1)-1)=1/1.1;
    Xi_hat_mnjf=ones(FI_J,N);
    
    elseif qq == 2 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 2:
    %
    % A +10% productivity shock to all firms to all countries except France.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'WORLDprod_p10'
    
    a_hat = ones(J*N,1);
    a_hat=a_hat/1.1;
    a_hat(start(15):start(15+1)-1)=1;   % Resets France to no shock
    Xi_hat_mnjf=ones(FI_J,N);
    
    elseif qq == 3
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 3:
    %
    % A +10% preference demand to all firms from Germany. Germany is the 10th
    % country in the WIOD.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'DEUpref_p10'
    
     Xi_hat_mnjf=ones(FI_J,N);
    Xi_hat_mnjf(:,10) = 1.1^(str2num(rhoT)-1);
    a_hat = ones(J*N,1);
    
    else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 4:
    %
    % A +10% preference demand to all firms from the World. France is the 15th
    % country in the WIOD.
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = 'WORLDpref_p10'
    
    Xi_hat_mnjf = ones(FI_J,N)*1.1^(str2num(rhoT)-1);
    Xi_hat_mnjf(:,15) = 1;
    a_hat = ones(J*N,1);
    end
  
    %% Running the Algorithm
   Pi_hat_n=ones(1,N);
   
   disp('BEGINNING COUNTERFACTUAL EXERCISE')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    % Step_4_Algorithm: same algorithm as Step_4_Algorithm_baseline.m, but now
    % we are able to feed in the model the calculated t+1 varialbes as t,
    % starting at the new equilibrium with no wedges
    % Load in the data after completing the convergence the first time, 
    % therefore calculating the new hat variables when the wedge is equal to 0. 
    % To do counterfactual exercises we start from this point, therefore
    % allowing our variables at time t to be equal to the hat variables
    % calculated previously. 
        alphaT = 4;
rhoT = '3';
J = 32;
rho = 3*ones(J,1);
sigma = 1.5*ones(J,1);
sigmaT = '1.5';
lambdaT = '1.000001';
lambda = 1.000001;
etaT = '1.000001';
eta = 1.000001;
varphi = 3;
  
       Step_4_Algorithm_counterfactual_CES_approx_Profit
% %     
% %     % Save the counterfactual output so it can be used for post analysis
% % 
      outputtable = [' X_mnj X_hat_mnj pi_mnjf1 pi_mnj1 ' ...
          'pi_c_mnj1 pi_c_nj1 pi_M1 pi_l1 pi_M_f1 pi_l_f1 '...
          'PC_hat_n1 PC_n P_hat_n P_hat_nj P_hat_mnj w_hat ' ...
          'a_hat a_hat_f Xi_hat_mnj Xi_hat_mnjf Pi_hat_n '...
          'final_goods1 intermediates1'] ;
% %  
%  
     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
        '_rho',rhoT,'_',shockT,'_Profit_approx.mat',outputtable,' -v7.3;')]);   

   display('STEP 4 COUNTERFACTUAL COMPLETED SUCCESFULLY')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % Step_5_counterfactual_output: Allows the user to easily analyse the
    %      counterfactual exercises produced before. Runs the CounterOutput.m
    %      function after specifying alpha type (labor share), gamma
    %      calibration and the shock scheme to be analysed.

                %% Load data

     sigmaT = '1.5';       
            % Time 0 data
           eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),...
    '_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
     '_rho',rhoT,'_Profit_approx.mat')]);

             % Time 1 data

                         eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),...
     '_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
     '_rho',rhoT,'_',shockT,'_Profit_approx.mat')]);

  
     Step_5_counterfactual_output_CES_approx_Profit
    
    % Save the counterfactual output so it can be used for post analysis

      outputtable = [' VA_hat_n VA_hat_mj VA_hat_mnj VA_hat_fj VA_hat_fnj '...
          'RGDP_def_hat_n RGDP_def_hat_fj RGDP_cpi_hat_n RGDP_cpi_hat_fj '...
          'RGDP_doubledef_hat_n RGDP_doubledef_hat_fj '...
          'real_wage_cpi_hat real_wage_def_hat real_wage_doubledef_hat '...
          'real_wage_income_cpi_hat_n real_wage_income_def_hat_n real_wage_income_doubledef_hat_n '...
          'TFP_hat_def_n TFP_hat_cpi_n TFP_hat_doubledef_n TFP_hat_cpi_france_j '...
          'imports_over_GDP_all imports_over_GDP_france_sector '...
          'exports_over_GDP_all exports_over_GDP_france_sector '...
          'GDP_deflator_n DoubleDeflation_deflator_n P_hat_n VA_fj_0 VA_mj_0'] ;
      %'TFP_hat_def_france_j TFP_hat_doubledef_france_j '...
%  
    sigmaT = '1.5';
     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_Profit_approx.mat',outputtable,' -v7.3;')]); 
  
%     
     disp('STEP 5 COMPLETED SUCCESFULLY')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step_6_Elastisticity_output: perform covariance decomposition at firm and
    %        sector level

    % Load files: only needed if run separately from Step 5
     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_Profit_approx.mat',outputtable,' -v7.3;')]); 

    % Shock size
    
   delta = 0.1;
    
   Step_6_Elasticity_output_CES_Profit

       outputtable = [' e_Y_FRf_cpi m_e_f_cpi cov_e_f_cpi check_f_cpi '...
           'e_Y_FRs_cpi m_e_mjFRs_cpi cov_e_mjFRs_cpi check_mjFRs_cpi '...
           'e_Y_s_cpi m_e_mj_cpi cov_e_mj_cpi check_mj_cpi '...
           'e_Y_FRf_def m_e_f_def cov_e_f_def check_f_def '...
           'e_Y_FRs_def m_e_mjFRs_def cov_e_mjFRs_def check_mjFRs_def '...
           'e_Y_s_def m_e_mj_def cov_e_mj_def check_mj_def '...
           'e_Y_FRf_doubledef m_e_f_doubledef cov_e_f_doubledef check_f_doubledef '...
           'e_Y_FRs_doubledef m_e_mjFRs_doubledef cov_e_mjFRs_doubledef check_mjFRs_doubledef '...
           'e_Y_s_doubledef m_e_mj_doubledef cov_e_mj_doubledef check_mj_doubledef'] ;
%  
     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_6_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_Profit_approx.mat',outputtable,' -v7.3;')]); 
     display('STEP 6 COMPLETED SUCCESFULLY')
    
    disp('COUNTERFACTUAL EXERCISE COMPLETED SUCCESFULLY')
end